This last test (ie: test_neofox) writes its output to a file named neofox/tests/resources/output_yyyymmddHHMMSS.txt. If there is an existing file named neofox/tests/resources/output_previous.txt then it loads both files in memory and compares them. It outputs whether there are some lost or gained columns and for the common columns it evaluates if the values are the same. If they are the same the file output_previous.txt is overwritten by the new file, otherwise it outputs the details of the differing columns.
Logs are written to the standard error and to the output folder by default. Optionally they can be written to a file by setting the environment variable NEOFOX_LOGFILE pointing to the desired file.
tar -xvf netMHCpan-4.0a.Linux.tar.gz
cd netMHCpan-4.0
wget http://www.cbs.dtu.dk/services/NetMHCpan-4.0/data.Linux.tar.gz
tar -xvf data.Linux.tar.gz
cd ..
export NEOFOX_NETMHCPAN=`pwd`/netMHCpan-4.0/netMHCpan
Configure NetMHCpan as explained in the file netMHCpan-4.0/netMHCpan-4.0.readme
tar -xvf netMHCIIpan-3.2.Linux.tar.gz
cd netMHCIIpan-3.2
# download the data
wget http://www.cbs.dtu.dk/services/NetMHCIIpan-3.2/data.Linux.tar.gz
tar -xvf data.Linux.tar.gz
cd ..
export NEOFOX_NETMHC2PAN=`pwd`/netMHCIIpan-3.2/netMHCIIpan
# install tcsh shell interpreter if not available yet
sudo apt-get install tcsh
Configure NetMHCpan as explained in the file netMHCIIpan-3.2/netMHCIIpan-3.2.readme
Unless indicated to the installer by flag --install-r-dependencies you will need to install manually some R dependencies. These dependencies are the following:
# Generated by the protocol buffer compiler. DO NOT EDIT!# sources: neoantigen.proto# plugin: python-betterprotofromdataclassesimportdataclassfromtypingimportListimportbetterproto
[docs]classZygosity(betterproto.Enum):"""*The zygosity of a given gene"""# *Two equal copies of the geneHOMOZYGOUS=0# *Two different copies of the geneHETEROZYGOUS=1# *Only one copy of the geneHEMIZYGOUS=2# *No copy of the geneLOSS=3
[docs]classMhc1Name(betterproto.Enum):"""*Valid names for MHC I classic genes"""A=0B=1C=2
[docs]classMhc2GeneName(betterproto.Enum):""" *Valid names for MHC II classic genes.DRA is not included in this list as it does not have much variability in the population and for our purpose isconsidered constant. """DRB1=0DPA1=1DPB1=2DQA1=3DQB1=4
[docs]classMhc2Name(betterproto.Enum):"""*Valid names for MHC II classic molecules"""DR=0DP=1DQ=2
[docs]@dataclassclassTranscript(betterproto.Message):# *The transcript identifier to which this neoepitope definition refers# (e.g.: Ensembl transcript id)identifier:str=betterproto.string_field(1)# *The genome assembly to which the gene definition refers to (e.g.: GRCh37,# GRCh38)assembly:str=betterproto.string_field(2)# *The gene symbol or gene identifier, optional as the transcript# unequivocally identifies a genegene:str=betterproto.string_field(3)
[docs]@dataclassclassMutation(betterproto.Message):# *The aminoacid position within the protein. 1-based, starting in the# N-terminusposition:int=betterproto.int32_field(1)# *Sequence of aminoacids for the wild type xmerwild_type_xmer:str=betterproto.string_field(2)# *IUPAC code for the wild type aminoacid in this positionwild_type_aminoacid:str=betterproto.string_field(3)# *Sequence of aminoacids for the mutated xmermutated_xmer:str=betterproto.string_field(4)# *IUPAC code for the mutated aminoacid in this positionmutated_aminoacid:str=betterproto.string_field(5)# *The left flanking region of variable size in IUPAC codesleft_flanking_region:str=betterproto.string_field(6)# *The size of the left flanking regionsize_left_flanking_region:int=betterproto.int32_field(7)# *The right flanking region of variable size in IUPAC codesright_flanking_region:str=betterproto.string_field(8)# *The size of the right flanking regionsize_right_flanking_region:int=betterproto.int32_field(9)
[docs]@dataclassclassNeoantigen(betterproto.Message):"""*A neoantigen minimal definition"""# *A unique identifier of a neoantigenidentifier:str=betterproto.string_field(1)# *Patient identifierpatient_identifier:str=betterproto.string_field(2)# *The transcript where the neoepitope was observedtranscript:"Transcript"=betterproto.message_field(3)# *The mutationmutation:"Mutation"=betterproto.message_field(4)# *Clonality estimation. At the moment this is a boolean indicating whether# there is clonality or not,there is no quantitive measurement at the moment.clonality_estimation:bool=betterproto.bool_field(5)# *Expression value of the transcript from RNA data. Range [0, +inf].rna_expression:float=betterproto.float_field(6)# *Variant allele frequency from the DNA. Range [0.0, 1.0]dna_variant_allele_frequency:float=betterproto.float_field(7)# *Variant allele frequency from the RNA. Range [0.0, 1.0]rna_variant_allele_frequency:float=betterproto.float_field(8)
[docs]@dataclassclassPatient(betterproto.Message):""" *The metadata required for analysis for a given patient + its patient identifier """# *Patient identifieridentifier:str=betterproto.string_field(1)# *Is RNA expression available?is_rna_available:bool=betterproto.bool_field(2)# *MHC I classic moleculesmhc1:List["Mhc1"]=betterproto.message_field(3)# *MHC II classic moleculesmhc2:List["Mhc2"]=betterproto.message_field(4)
[docs]@dataclassclassAnnotation(betterproto.Message):"""*This is a generic class to hold annotations from Neofox"""# *The name of the annotationname:str=betterproto.string_field(1)# *The value of the annotationTODO: make this better to hold int, float and# boolvalue:str=betterproto.string_field(2)
[docs]@dataclassclassNeoantigenAnnotations(betterproto.Message):"""*A set of annotations for a neoantigen"""# *A unique identifier of a neoantigenneoantigen_identifier:str=betterproto.string_field(1)# *List of annotationsannotations:List["Annotation"]=betterproto.message_field(2)# *The annotatorannotator:str=betterproto.string_field(3)# *The version of the annotatorannotator_version:str=betterproto.string_field(4)# *A timestamp determined when the annotation was createdtimestamp:str=betterproto.string_field(5)# *Annotation resources MD5 hashresources_hash:str=betterproto.string_field(6)
[docs]@dataclassclassMhc1(betterproto.Message):""" *Models MHC I alleles related to the same MHC I gene, i.e. 2 alleles/2 isoforms per gene """# *MHC I gene namename:"Mhc1Name"=betterproto.enum_field(1)# *Zygosity of the genezygosity:"Zygosity"=betterproto.enum_field(2)# *The alleles of the gene (0, 1 or 2)alleles:List["MhcAllele"]=betterproto.message_field(3)
[docs]@dataclassclassMhc2(betterproto.Message):""" *Models MHC II alleles related to the same MHC II protein, i.e. 4 isoforms related to 2 genes with 2 alleles each """# *MHC II molecule namename:"Mhc2Name"=betterproto.enum_field(1)# *List of MHC II genesgenes:List["Mhc2Gene"]=betterproto.message_field(2)# *Different combinations of MHC II alleles building different isoformsisoforms:List["Mhc2Isoform"]=betterproto.message_field(3)
[docs]@dataclassclassMhc2Isoform(betterproto.Message):"""*MHC II isoform"""# *Name to refer to the MHC II isoformname:str=betterproto.string_field(1)# *The alpha chain of the isoformalpha_chain:"MhcAllele"=betterproto.message_field(2)# *The beta chain of the isoformbeta_chain:"MhcAllele"=betterproto.message_field(3)
[docs]@dataclassclassMhc2Gene(betterproto.Message):"""*MHC II gene"""# *MHC II gene namename:"Mhc2GeneName"=betterproto.enum_field(1)# *Zygosity of the genezygosity:"Zygosity"=betterproto.enum_field(2)# *The alleles of the gene (0, 1 or 2)alleles:List["MhcAllele"]=betterproto.message_field(3)
[docs]@dataclassclassMhcAllele(betterproto.Message):""" *MHC allele representation. It does not include non synonymous changes to the sequence, changes in the non coding regionor changes in expression. See http://hla.alleles.org/nomenclature/naming.html for details """# *HLA full name as provided by the user (e.g.: HLA-DRB1*13:01:02:03N). This# will be parsed into name, gene and group.Any digit format is allowed for# this field (ie: 4, 6 or 8 digits), 2 digits names are not specific enough# for ourpurpose and thus invalidfull_name:str=betterproto.string_field(1)# *A specific HLA protein (e.g. HLA-DRB1*13:01). Alleles whose numbers differ# in group and protein must differ in oneor more nucleotide substitutions# that change the amino acid sequence of the encoded protein.This name is# normalized to avoid different representations of the same allele. For# instance both HLA-DRB113:01 andHLA-DRB1*13:01:02:03N will be transformed# into their normalised version HLA-DRB1*13:01. This name is also truncatedto# 4 digits. 2 digits names are not specific enough for our purpose and thus# invalidname:str=betterproto.string_field(2)# *The gene from either MHC I or II (e.g. DRB1, A) (this information is# redundant with the Mhc1Gene.name andMhc2Gene.name but it is convenient to# have this at this level too, code will check for data coherence)gene:str=betterproto.string_field(3)# *A group of alleles defined by a common serotype ie: Serological antigen# carried by an allotype (e.g. 13 from HLA-DRB1*13)group:str=betterproto.string_field(4)# *A specific protein (e.g.: 02 from HLA-DRB1*13:02)protein:str=betterproto.string_field(5)
## Copyright (c) 2020-2030 Translational Oncology at the Medical Center of the Johannes Gutenberg-University Mainz gGmbH.## This file is part of Neofox# (see https://github.com/tron-bioinformatics/neofox).## This program is free software: you can redistribute it and/or modify# it under the terms of the GNU General Public License as published by# the Free Software Foundation, either version 3 of the License, or# (at your option) any later version.## This program is distributed in the hope that it will be useful,# but WITHOUT ANY WARRANTY; without even the implied warranty of# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the# GNU General Public License for more details.## You should have received a copy of the GNU General Public License# along with this program. If not, see <http://www.gnu.org/licenses/>.#fromtypingimportListimportpandasaspdimportbetterprotofrombetterprotoimportCasingfrompandas.io.jsonimportjson_normalizeimportreimportdifflibfromcollectionsimportdefaultdictimportjsonfromneofox.model.neoantigenimportNeoantigen,Transcript,Mutation,Patient,NeoantigenAnnotations,Mhc2Name, \
Mhc2GeneName,Zygosity,Mhc2Gene,Mhc2,Mhc2Isoform,MhcAllele,Mhc1Name,Mhc1FIELD_SUBSTITUTION='substitution'FIELD_VAF_DNA='VAF_in_tumor'FIELD_VAF_RNA='VAF_in_RNA'FIELD_RNA_EXPRESSION='VAF_RNA_raw'FIELD_TRANSCRIPT_EXPRESSION='transcript_expression'FIELD_TRANSCRIPT='UCSC_transcript'FIELD_GENE='gene'FIELD_WILD_TYPE_XMER='[WT]_+-13_AA_(SNV)_/_-15_AA_to_STOP_(INDEL)'FIELD_MUTATED_XMER='+-13_AA_(SNV)_/_-15_AA_to_STOP_(INDEL)'
[docs]@staticmethoddefparse_candidate_file(candidate_file:str,patient_id:str=None)->List[Neoantigen]:""" :param candidate_file: the path to an neoantigen candidate input file :param patient_id: the patient identifier for all neoantigens in the input file, if not provided it is expected as column named `patient.id` or `patient` :return neoantigens in model objects """data=pd.read_csv(candidate_file,sep='\t')# filter out indelsdata=data[~data['substitution'].isna()]ModelConverter._enrich_candidate_table(data)neoantigens=[]for_,candidate_entryindata.iterrows():neoantigens.append(ModelConverter._candidate_entry2model(candidate_entry,patient_id=patient_id))returnneoantigens
[docs]@staticmethoddefparse_patients_file(patients_file:str)->List[Patient]:""" :param patients_file: the file to patients data CSV file :return: the parsed CSV into model objects """split_comma_separated_list=lambdax:x.split(',')df=pd.read_csv(patients_file,sep='\t',converters={'mhcIAlleles':split_comma_separated_list,'mhcIIAlleles':split_comma_separated_list})returnModelConverter.patient_metadata_csv2objects(df)
[docs]@staticmethoddefparse_neoantigens_file(neoantigens_file:str)->List[Neoantigen]:""" :param neoantigens_file: the file to neoantigens data CSV file :return: the parsed CSV into model objects """returnModelConverter.neoantigens_csv2objects(pd.read_csv(neoantigens_file,sep='\t').fillna(""))
[docs]@staticmethoddefobjects2dataframe(model_objects:List[betterproto.Message])->pd.DataFrame:""" :param model_objects: list of objects of subclass of betterproto.Message """returnjson_normalize(data=[n.to_dict(include_default_values=True)forninmodel_objects])
[docs]@staticmethoddefobjects2json(model_objects:List[betterproto.Message],output_file:str):""" :param model_objects: list of objects of subclass of betterproto.Message """withopen(output_file,"w")asf:json.dump([o.to_dict(casing=Casing.SNAKE)foroinmodel_objects],f)
[docs]@staticmethoddefobject2series(model_object:betterproto.Message)->pd.Series:""" :param model_object: object of subclass of betterproto.Message """returnjson_normalize(data=model_object.to_dict(casing=Casing.SNAKE,include_default_values=True)).iloc[0]
[docs]@staticmethoddefobject2flat_dict(model:betterproto.Message)->dict:""" Transforms a model object into a flat dict. Nested fields are concatenated with a dot """returnModelConverter.object2series(model).to_dict()
[docs]@staticmethoddefneoantigens_csv2object(series:pd.Series)->Neoantigen:"""transforms an entry from a CSV into an object"""returnNeoantigen().from_dict(ModelConverter._flat_dict2nested_dict(flat_dict=series.to_dict()))
[docs]@staticmethoddefpatient_metadata_csv2objects(dataframe:pd.DataFrame)->List[Patient]:"""transforms an patients CSV into a list of objects"""patients=[]for_,rowindataframe.iterrows():patient_dict=row.to_dict()patient=Patient().from_dict(patient_dict)patient.mhc1=ModelConverter.parse_mhc1_alleles(patient_dict['mhcIAlleles'])patient.mhc2=ModelConverter.parse_mhc2_alleles(patient_dict['mhcIIAlleles'])patients.append(patient)returnpatients
[docs]@staticmethoddefneoantigens_csv2objects(dataframe:pd.DataFrame)->List[Neoantigen]:"""transforms an patients CSV into a list of objects"""neoantigens=[]for_,rowindataframe.iterrows():neoantigens.append(Neoantigen().from_dict(ModelConverter._flat_dict2nested_dict(flat_dict=row.to_dict())))returnneoantigens
[docs]@staticmethoddefannotations2tall_skinny_table(neoantigen_annotations:List[NeoantigenAnnotations])->pd.DataFrame:dfs=[]fornainneoantigen_annotations:df=pd.DataFrame([a.to_dict()forainna.annotations])df['neoantigen_identifier']=na.neoantigen_identifierdfs.append(df[df.value!="NA"])# avoid writing NA valuesreturnpd.concat(dfs)
@staticmethoddef_flat_dict2nested_dict(flat_dict:dict)->dict:"""transforms a flattened dict into a nested dict, assuming that the dot indicates a nested level"""nested_dict=defaultdict(lambda:{})fork,vinflat_dict.items():splitted_k=k.split('.')iflen(splitted_k)>2:raiseNotImplemented("Support for dictionaries nested more than one level is not implemented")iflen(splitted_k)==2:nested_dict[splitted_k[0]][splitted_k[1]]=velse:nested_dict[k]=vreturndict(nested_dict)@staticmethoddef_candidate_entry2model(candidate_entry:dict,patient_id:str)->Neoantigen:"""parses an row from a candidate file into a model object"""transcript=Transcript()transcript.assembly='hg19'transcript.gene=candidate_entry.get(FIELD_GENE)transcript.identifier=candidate_entry.get(FIELD_TRANSCRIPT)mutation=Mutation()mutation.position=candidate_entry.get('position')mutation.wild_type_aminoacid=candidate_entry.get('wild_type_aminoacid')mutation.mutated_aminoacid=candidate_entry.get('mutated_aminoacid')mutation.left_flanking_region=candidate_entry.get('left_flanking_region')mutation.right_flanking_region=candidate_entry.get('right_flanking_region')neoantigen=Neoantigen()neoantigen.patient_identifier=patient_idifpatient_idelsecandidate_entry.get('patient')neoantigen.mutation=mutationneoantigen.transcript=transcript# clonality estimation is not present in candidate file at the momentneoantigen.clonality_estimation=None# missing RNA expression values are represented as -1vaf_rna_raw=candidate_entry.get(FIELD_TRANSCRIPT_EXPRESSION)neoantigen.rna_expression=vaf_rna_rawifvaf_rna_raw>=0elseNonevaf_in_rna=candidate_entry.get(FIELD_VAF_RNA)neoantigen.rna_variant_allele_frequency=vaf_in_rnaifvaf_in_rna>=0elseNoneneoantigen.dna_variant_allele_frequency=candidate_entry.get(FIELD_VAF_DNA)returnneoantigen@staticmethoddef_enrich_candidate_table(data:pd.DataFrame):"""parses some data from the candidate table into the right fields in the model objects"""data['wild_type_aminoacid']=data[FIELD_SUBSTITUTION].transform(lambdax:re.search("(\w)\d+\w",x).group(1))data['mutated_aminoacid']=data[FIELD_SUBSTITUTION].transform(lambdax:re.search("\w\d+(\w)",x).group(1))data['position']=data[FIELD_SUBSTITUTION].transform(lambdax:int(re.search("\w(\d+)\w",x).group(1)))data['left_flanking_region']=data[[FIELD_MUTATED_XMER,FIELD_WILD_TYPE_XMER]].apply(lambdax:ModelConverter._get_matching_region(x[0],x[1]),axis=1)data['right_flanking_region']=data[[FIELD_MUTATED_XMER,FIELD_WILD_TYPE_XMER]].apply(lambdax:ModelConverter._get_matching_region(x[0],x[1],match=1),axis=1)@staticmethoddef_get_matching_region(sequence1:str,sequence2:str,match:int=0)->str:"""fetches an aligned block within two sequences"""match=difflib.SequenceMatcher(None,sequence1,sequence2).get_matching_blocks()[match]returnsequence1[match.a:match.a+match.size]
[docs]@staticmethoddefparse_mhc1_alleles(alleles:List[str])->List[Mhc1]:isoforms=[]parsed_alleles=list(map(ModelConverter.parse_mhc_allele,alleles))ModelConverter._validate_mhc1_alleles(parsed_alleles)# do we need to validate genes anymore? add test creating MhcAllele with bad gene and see what happensforgene_nameinMhc1Name:gene_alleles=list(filter(lambdaa:a.gene==gene_name.name,parsed_alleles))zygosity=ModelConverter._get_zygosity_from_alleles(gene_alleles)ifzygosity==Zygosity.HOMOZYGOUS:gene_alleles=[gene_alleles[0]]# we don't want repeated instances of the same alleleisoforms.append(Mhc1(name=gene_name,zygosity=zygosity,alleles=gene_alleles))returnisoforms
[docs]@staticmethoddefparse_mhc2_alleles(alleles:List[str])->List[Mhc2]:mhc2s=[]parsed_alleles=list(map(ModelConverter.parse_mhc_allele,alleles))ModelConverter._validate_mhc2_alleles(parsed_alleles)# do we need to validate genes anymore? add test creating MhcAllele with bad gene and see what happensforisoform_nameinMhc2Name:isoform_alleles=list(filter(lambdaa:isoform_name.nameina.gene,parsed_alleles))genes=[]forgene_nameinModelConverter.GENES_BY_MOLECULE.get(isoform_name):gene_alleles=list(filter(lambdaa:a.gene==gene_name.name,isoform_alleles))zygosity=ModelConverter._get_zygosity_from_alleles(gene_alleles)ifzygosity==Zygosity.HOMOZYGOUS:gene_alleles=[gene_alleles[0]]# we don't want repeated instances of the same allelegenes.append(Mhc2Gene(name=gene_name,zygosity=zygosity,alleles=gene_alleles))isoforms=ModelConverter._get_mhc2_isoforms(isoform_name,genes)mhc2s.append(Mhc2(name=isoform_name,genes=genes,isoforms=isoforms))returnmhc2s
@staticmethoddef_get_mhc2_isoforms(isoform_name:Mhc2Name,genes:List[Mhc2Gene])->List[Mhc2Isoform]:isoforms=[]ifisoform_name==Mhc2Name.DR:assertlen(genes)<=1,"More than one gene provided for MHC II DR"# alpha chain of the MHC II DR is not modelled as it is constantisoforms=[Mhc2Isoform(name=a.name,alpha_chain=None,beta_chain=a)forgingenesforaing.alleles]elifisoform_name==Mhc2Name.DP:assertlen(genes)<=2,"More than two genes provided for MHC II DP"alpha_alleles=[aforgingenesifg.name==Mhc2GeneName.DPA1foraing.alleles]beta_alleles=[aforgingenesifg.name==Mhc2GeneName.DPB1foraing.alleles]isoforms=[Mhc2Isoform(name=ModelConverter.get_mhc2_isoform_name(a,b),alpha_chain=a,beta_chain=b)forainalpha_allelesforbinbeta_alleles]elifisoform_name==Mhc2Name.DQ:assertlen(genes)<=2,"More than two genes provided for MHC II DQ"alpha_alleles=[aforgingenesifg.name==Mhc2GeneName.DQA1foraing.alleles]beta_alleles=[aforgingenesifg.name==Mhc2GeneName.DQB1foraing.alleles]isoforms=[Mhc2Isoform(name=ModelConverter.get_mhc2_isoform_name(a,b),alpha_chain=a,beta_chain=b)forainalpha_allelesforbinbeta_alleles]returnisoforms
[docs]@staticmethoddefget_mhc2_isoform_name(a:MhcAllele,b:MhcAllele):# NOTE: this is needed as jus setting alpha chain to None wouldn't work with protobufifaisnotNoneanda.name:return"{}-{}".format(a.name,b.name.replace("HLA-",""))else:returnb.name
[docs]@staticmethoddefparse_mhc_allele(allele:str)->MhcAllele:# infers gene, group and protein from the namematch=ModelConverter.HLA_ALLELE_PATTERN.match(allele)assertmatchisnotNone,"Allele does not match HLA allele pattern {}".format(allele)gene=match.group(1)group=match.group(2)protein=match.group(3)# builds a normalized representation of the allelename="HLA-{gene}*{serotype}:{protein}".format(gene=gene,serotype=group,protein=protein)# ensures that full name stores the complete allele as provided but normalizes# its representationfull_name=namesix_digits_id=match.group(4)ifsix_digits_idisnotNoneandsix_digits_id!="":full_name=full_name+":{}".format(six_digits_id)eight_digits_id=match.group(5)ifeight_digits_idisnotNoneandeight_digits_id!="":full_name=full_name+":{}".format(eight_digits_id)expression_change=match.group(6)ifexpression_changeisnotNoneandexpression_change!="":full_name=full_name+expression_changereturnMhcAllele(full_name=full_name,name=name,gene=gene,group=group,protein=protein)
@staticmethoddef_get_zygosity_from_alleles(alleles:List[MhcAllele])->Zygosity:assertlen(set([a.geneforainalleles]))<=1,"Trying to get zygosity from alleles of different genes"assertlen(alleles)<=2,"More than 2 alleles for gene {}".format(alleles[0].gene)iflen(alleles)==2andalleles[0].name==alleles[1].name:zygosity=Zygosity.HOMOZYGOUSeliflen(alleles)==2:zygosity=Zygosity.HETEROZYGOUSeliflen(alleles)==1:zygosity=Zygosity.HEMIZYGOUSelse:zygosity=Zygosity.LOSSreturnzygosity@staticmethoddef_validate_mhc1_alleles(parsed_alleles:List[MhcAllele]):forainparsed_alleles:asserta.geneinMhc1Name.__members__,"Gene from MHC I allele is not valid: {}".format(a.gene)@staticmethoddef_validate_mhc2_alleles(parsed_alleles:List[MhcAllele]):forainparsed_alleles:asserta.geneinMhc2GeneName.__members__,"Gene from MHC II allele is not valid: {}".format(a.gene)
## Copyright (c) 2020-2030 Translational Oncology at the Medical Center of the Johannes Gutenberg-University Mainz gGmbH.## This file is part of Neofox# (see https://github.com/tron-bioinformatics/neofox).## This program is free software: you can redistribute it and/or modify# it under the terms of the GNU General Public License as published by# the Free Software Foundation, either version 3 of the License, or# (at your option) any later version.## This program is distributed in the hope that it will be useful,# but WITHOUT ANY WARRANTY; without even the implied warranty of# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the# GNU General Public License for more details.## You should have received a copy of the GNU General Public License# along with this program. If not, see <http://www.gnu.org/licenses/>.#fromtypingimportListfromneofox.model.neoantigenimportAnnotation,Mhc2,Mhc2GeneName,MhcAlleleNOT_AVAILABLE_VALUE="NA"
[docs]@staticmethoddefbuild_annotation(name,value):ifisinstance(value,bool):# prints booleans as 0/1 stringsvalue="1"ifvalueelse"0"ifnotisinstance(value,str)andnotisinstance(value,type(None)):ifisinstance(value,float):value="{0:.5g}".format(round(value,5))else:value=str(value)ifvalueisNone:value=NOT_AVAILABLE_VALUEreturnAnnotation(name=name,value=value)
## Copyright (c) 2020-2030 Translational Oncology at the Medical Center of the Johannes Gutenberg-University Mainz gGmbH.## This file is part of Neofox# (see https://github.com/tron-bioinformatics/neofox).## This program is free software: you can redistribute it and/or modify# it under the terms of the GNU General Public License as published by# the Free Software Foundation, either version 3 of the License, or# (at your option) any later version.## This program is distributed in the hope that it will be useful,# but WITHOUT ANY WARRANTY; without even the implied warranty of# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the# GNU General Public License for more details.## You should have received a copy of the GNU General Public License# along with this program. If not, see <http://www.gnu.org/licenses/>.#importbase64importhashlibimportbetterprotofromneofox.model.conversionimportModelConverterfromneofox.exceptionsimportNeofoxDataValidationExceptionfromneofox.model.neoantigenimportNeoantigen,Mutation,Transcript,Patient,MhcAllele,Mhc2Isoform,Mhc1, \
Mhc1Name,Zygosity,Mhc2,Mhc2Name,Mhc2GeneNamefromBio.Alphabet.IUPACimportExtendedIUPACProtein,IUPACData
[docs]@staticmethoddefvalidate(model:betterproto.Message):# TODO: make this method capture appropriately validation issues whend dealing with int and floattry:model.__bytes__()exceptExceptionase:raiseNeofoxDataValidationException(e)
[docs]@staticmethoddefvalidate_neoantigen(neoantigen:Neoantigen)->Neoantigen:# checks format consistency firstModelValidator.validate(neoantigen)try:# checks gene# TODO: do we want to verify existence of gene and transcript id?neoantigen.transcript=ModelValidator._validate_transcript(neoantigen.transcript)# checks mutationneoantigen.mutation=ModelValidator._validate_mutation(neoantigen.mutation)# check the expression valuesModelValidator._validate_expression_values(neoantigen)exceptAssertionErrorase:raiseNeofoxDataValidationException(e)# infer other fields from the modelreturnModelValidator._enrich_neoantigen(neoantigen)
[docs]@staticmethoddefvalidate_patient(patient:Patient)->Patient:# checks format consistency firstModelValidator.validate(patient)try:# checks that patient id is not empty considering white spacespatient_id=patient.identifier.strip()ifpatient.identifierelsepatient.identifierassertpatient_idisnotNoneandpatient_id!="","Patient identifier is empty"patient.identifier=patient_id# TODO: validate new model with isoforms, genes and alleles# checks MHC Ivalidated_mhc1s=[]forminpatient.mhc1:validated_mhc1s.append(ModelValidator._validate_mhc1(m))patient.mhc1=validated_mhc1s# checks MHC IIvalidated_mhc2s=[]forminpatient.mhc2:validated_mhc2s.append(ModelValidator._validate_mhc2(m))patient.mhc2=validated_mhc2sexceptAssertionErrorase:raiseNeofoxDataValidationException(e)returnpatient
@staticmethoddef_validate_mhc1(mhc1:Mhc1)->Mhc1:assertmhc1.nameinMhc1Name,"Invalid MHC I name"assertmhc1.zygosityinZygosity,"Invalid zygosity"alleles=mhc1.allelesifmhc1.zygosityin[Zygosity.HOMOZYGOUS,Zygosity.HEMIZYGOUS]:assertlen(alleles)==1,"An homozygous gene must have 1 allele and not {}".format(len(alleles))elifmhc1.zygosity==Zygosity.HETEROZYGOUS:assertlen(alleles)==2,"An heterozygous or hemizygous gene must have 2 alleles and not {}".format(len(alleles))elifmhc1.zygosity==Zygosity.LOSS:assertlen(alleles)==0,"An lost gene must have 0 alleles and not {}".format(len(alleles))validated_alleles=[]foralleleinalleles:validated_allele=ModelValidator.validate_mhc_allele_representation(allele)validated_alleles.append(validated_allele)assertvalidated_allele.gene==mhc1.name.name, \
"The allele referring to gene {} is inside gene {}".format(validated_allele.gene,mhc1.name.name)mhc1.alleles=validated_allelesreturnmhc1@staticmethoddef_validate_mhc2(mhc2:Mhc2)->Mhc2:assertmhc2.nameinMhc2Name,"Invalid MHC II name"genes=mhc2.genesforgeneingenes:assertgene.nameinMhc2GeneName,"Invalid gene name from MHC II"assertgene.nameinModelConverter.GENES_BY_MOLECULE.get(mhc2.name), \
"Gene {} referring to isoform {}".format(gene.name,mhc2.name.name)assertgene.zygosityinZygosity,"Invalid zygosity"alleles=gene.allelesifgene.zygosity==Zygosity.HOMOZYGOUS:assertlen(alleles)==1,"An homozygous gene must have 1 allele and not {}".format(len(alleles))elifgene.zygosityin[Zygosity.HETEROZYGOUS,Zygosity.HEMIZYGOUS]:assertlen(alleles)==2,"An heterozygous or hemizygous gene must have 2 alleles and not {}".format(len(alleles))elifgene.zygosity==Zygosity.LOSS:assertlen(alleles)==0,"An lost gene must have 0 alleles and not {}".format(len(alleles))validated_alleles=[]foralleleinalleles:validated_allele=ModelValidator.validate_mhc_allele_representation(allele)validated_alleles.append(validated_allele)assertvalidated_allele.gene==gene.name.name, \
"The allele referring to gene {} is inside gene {}".format(validated_allele.gene,gene.name.name)gene.alleles=validated_allelesisoforms=mhc2.isoformsvalidated_isoforms=[]forisoforminisoforms:validated_isoform=ModelValidator.validate_mhc2_isoform_representation(isoform)validated_isoforms.append(validated_isoform)ifmhc2.name!=Mhc2Name.DR:assertvalidated_isoform.alpha_chain.namein[a.nameforgingenesforaing.alleles], \
"Alpha chain allele not present in th list of alleles"assertvalidated_isoform.beta_chain.namein[a.nameforgingenesforaing.alleles], \
"Beta chain allele not present in th list of alleles"mhc2.isoforms=validated_isoformsreturnmhc2
[docs]@staticmethoddefvalidate_mhc_allele_representation(allele:MhcAllele)->MhcAllele:try:full_name=Noneifallele.full_name:# infers gene, group and protein from the namematch=ModelConverter.HLA_ALLELE_PATTERN.match(allele.full_name)assertmatchisnotNone,"Allele does not match HLA allele pattern {}".format(allele.name)gene=match.group(1)group=match.group(2)protein=match.group(3)full_name=allele.full_nameelifallele.name:# infers gene, group and protein from the namematch=ModelConverter.HLA_ALLELE_PATTERN.match(allele.name)assertmatchisnotNone,"Allele does not match HLA allele pattern {}".format(allele.name)gene=match.group(1)group=match.group(2)protein=match.group(3)elifallele.geneandallele.groupandallele.protein:# infers name from gene, group and proteingene=allele.genegroup=allele.groupprotein=allele.proteinelse:raiseNeofoxDataValidationException("HLA allele missing required fields, either name or gene, group and protein must be provided")assertgeneinlist(Mhc1Name.__members__.keys())+list(Mhc2GeneName.__members__.keys()), \
"Gene not from classic MHC: {}".format(gene)# builds the final allele representation and validates it just in casename="HLA-{gene}*{serotype}:{protein}".format(gene=gene,serotype=group,protein=protein)match=ModelConverter.HLA_ALLELE_PATTERN.match(name)assertmatchisnotNone,"Allele does not match HLA allele pattern {}".format(name)exceptAssertionErrorase:raiseNeofoxDataValidationException(e)returnMhcAllele(full_name=full_nameiffull_nameelsename,name=name,gene=gene,group=group,protein=protein)
[docs]@staticmethoddefvalidate_mhc2_isoform_representation(isoform:Mhc2Isoform)->Mhc2Isoform:try:ifisoform.name:# infers alpha and beta chainsmatch=ModelConverter.HLA_MOLECULE_PATTERN.match(isoform.name)ifmatch:alpha_chain=ModelValidator.validate_mhc_allele_representation(MhcAllele(name=match.group(1)))beta_chain=ModelValidator.validate_mhc_allele_representation(MhcAllele(name=match.group(2)))else:match=ModelConverter.HLA_DR_MOLECULE_PATTERN.match(isoform.name)assertmatchisnotNone,"Molecule does not match HLA isoform pattern {}".format(isoform.name)alpha_chain=MhcAllele()beta_chain=ModelValidator.validate_mhc_allele_representation(MhcAllele(name=match.group(1)))elifisoform.alpha_chainandisoform.beta_chain:# infers name from gene, group and proteinalpha_chain=ModelValidator.validate_mhc_allele_representation(isoform.alpha_chain)beta_chain=ModelValidator.validate_mhc_allele_representation(isoform.beta_chain)else:raiseNeofoxDataValidationException("HLA isoform missing required fields")# builds the final allele representation and validates it just in casename=ModelConverter.get_mhc2_isoform_name(alpha_chain,beta_chain)match=ModelConverter.HLA_MOLECULE_PATTERN.match(name)match2=ModelConverter.HLA_DR_MOLECULE_PATTERN.match(name)assertmatchisnotNoneormatch2isnotNone, \
"Molecule does not match HLA isoform pattern {}".format(name)exceptAssertionErrorase:raiseNeofoxDataValidationException(e)returnMhc2Isoform(name=name,alpha_chain=alpha_chain,beta_chain=beta_chain)
@staticmethoddef_validate_expression_values(neoantigen):assertneoantigen.rna_expressionisNoneorneoantigen.rna_expression>=0, \
"RNA expression should be a positive integer or zero {}".format(neoantigen.rna_expression)ModelValidator._validate_vaf(neoantigen.dna_variant_allele_frequency)ModelValidator._validate_vaf(neoantigen.rna_variant_allele_frequency)@staticmethoddef_validate_mutation(mutation:Mutation)->Mutation:# checks aminoacidsmutation.mutated_aminoacid=ModelValidator._validate_aminoacid(mutation.mutated_aminoacid)mutation.wild_type_aminoacid=ModelValidator._validate_aminoacid(mutation.wild_type_aminoacid)# checks left and right flanking regionsassertmutation.left_flanking_regionisnotNone,"Empty left flanking region"mutation.left_flanking_region=mutation.left_flanking_region.strip()assertlen(mutation.left_flanking_region)>0,"Empty left flanking region"foraainmutation.left_flanking_region:ModelValidator._validate_aminoacid(aa)assertmutation.right_flanking_regionisnotNone,"Empty right flanking region"mutation.right_flanking_region=mutation.right_flanking_region.strip()assertlen(mutation.right_flanking_region)>0,"Empty right flanking region"foraainmutation.right_flanking_region:ModelValidator._validate_aminoacid(aa)# checks the positionassertmutation.positionisnotNone,"Empty position"assertisinstance(mutation.position,int),"Position must be an integer"assertmutation.position>0,"Position must be a 1-based positive integer"returnmutation@staticmethoddef_validate_transcript(transcript:Transcript)->Transcript:# TODO: validate that gene symbol existsgene_name=transcript.gene.strip()iftranscript.geneelsetranscript.geneassertgene_nameisnotNoneandlen(gene_name)>0,"Empty gene symbol"transcript.gene=gene_name# TODO: validate that transcript identifier existstranscript_identifier=transcript.identifier.strip()iftranscript.identifierelsetranscript.identifierasserttranscript_identifierisnotNoneandlen(transcript_identifier)>0,"Empty transcript identifier"transcript.identifier=transcript_identifier# TODO: support other assembliesassembly=transcript.assemblyiftranscript.assemblyelse"hg19"assertassembly=="hg19","Other reference genome than hg19 is not supported"transcript.assembly=assemblyreturntranscript@staticmethoddef_validate_vaf(vaf):assertvafisNoneor0.0<=vaf<=1.0,"VAF should be a positive integer or zero {}".format(vaf)@staticmethoddef_enrich_neoantigen(neoantigen:Neoantigen)->Neoantigen:neoantigen.mutation.wild_type_xmer="".join([neoantigen.mutation.left_flanking_region,neoantigen.mutation.wild_type_aminoacid,neoantigen.mutation.right_flanking_region])neoantigen.mutation.mutated_xmer="".join([neoantigen.mutation.left_flanking_region,neoantigen.mutation.mutated_aminoacid,neoantigen.mutation.right_flanking_region])neoantigen.mutation.size_left_flanking_region=len(neoantigen.mutation.left_flanking_region)neoantigen.mutation.size_right_flanking_region=len(neoantigen.mutation.right_flanking_region)neoantigen.identifier=ModelValidator.generate_neoantigen_identifier(neoantigen)returnneoantigen@staticmethoddef_validate_aminoacid(aminoacid):assertaminoacidisnotNone,"Aminoacid field cannot be empty"aminoacid=aminoacid.strip()assertisinstance(aminoacid,str),"Aminoacid has to be a string"iflen(aminoacid)==3:assertaminoacidinIUPACData.protein_letters_3to1_extended.keys(), \
"Non existing 3 letter aminoacid {}".format(aminoacid)aminoacid=IUPACData.protein_letters_3to1_extended.get(aminoacid)iflen(aminoacid)==1:assertaminoacidinExtendedIUPACProtein.letters,"Non existing aminoacid {}".format(aminoacid)else:assertFalse,"Invalid aminoacid {}".format(aminoacid)returnaminoacid
[docs]@staticmethoddefgenerate_neoantigen_identifier(neoantigen:Neoantigen)->str:neoantigen.identifier=None# this needs to be done otherwise we cannot rreproduce the id after it is setreturnbase64.b64encode(hashlib.md5(neoantigen.to_json().encode('utf8')).digest()).decode('utf8')
#!/usr/bin/env python## Copyright (c) 2020-2030 Translational Oncology at the Medical Center of the Johannes Gutenberg-University Mainz gGmbH.## This file is part of Neofox# (see https://github.com/tron-bioinformatics/neofox).## This program is free software: you can redistribute it and/or modify# it under the terms of the GNU General Public License as published by# the Free Software Foundation, either version 3 of the License, or# (at your option) any later version.## This program is distributed in the hope that it will be useful,# but WITHOUT ANY WARRANTY; without even the implied warranty of# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the# GNU General Public License for more details.## You should have received a copy of the GNU General Public License# along with this program. If not, see <http://www.gnu.org/licenses/>.#importloggingimportosimporttimefromtypingimportListimportlogzerofromlogzeroimportloggerfromdask.distributedimportClientfromneofox.references.referencesimportReferenceFolder,AvailableAlleles,DependenciesConfigurationfromneofoximportNEOFOX_LOG_FILE_ENVfromneofox.annotatorimportNeoantigenAnnotatorfromneofox.exceptionsimportNeofoxConfigurationException,NeofoxDataValidationExceptionfromneofox.model.neoantigenimportNeoantigenAnnotations,Neoantigen,Patientfromneofox.model.validationimportModelValidator
[docs]classNeoFox:def__init__(self,neoantigens:List[Neoantigen],patient_id:str,patients:List[Patient],num_cpus:int,work_folder=None,output_prefix=None,reference_folder:ReferenceFolder=None,configuration:DependenciesConfiguration=None):# initialise logsself._initialise_logs(output_prefix,work_folder)# initialise dask# TODO: number of threads is hard coded. Is there a better value for this?self.dask_client=Client(processes=True,n_workers=num_cpus,threads_per_worker=4)# intialize references folder and configuration# NOTE: uses the reference folder and config passed as a parameter if exists, this is here to make it# testable with fake objectsself.reference_folder=reference_folderifreference_folderelseReferenceFolder()self.configuration=configurationifconfigurationelseDependenciesConfiguration()ifneoantigensisNoneorlen(neoantigens)==0orpatientsisNoneorlen(patients)==0:raiseNeofoxConfigurationException("Missing input data to run Neofox")# TODO: avoid overriding patient id parameterforninneoantigens:ifn.patient_identifierisNone:n.patient_identifier=patient_id# validates input dataself.neoantigens=[ModelValidator.validate_neoantigen(n)forninneoantigens]self.patients={patient.identifier:ModelValidator.validate_patient(patient)forpatientinpatients}self._validate_input_data()logger.info("Data loaded")def_initialise_logs(self,output_prefix,work_folder):ifwork_folderandos.path.exists(work_folder):logfile=os.path.join(work_folder,"{}.log".format(output_prefix))else:logfile=os.environ.get(NEOFOX_LOG_FILE_ENV)iflogfileisnotNone:logzero.logfile(logfile)# TODO: this does not worklogzero.loglevel(logging.DEBUG)logger.info("Loading data...")def_validate_input_data(self):patient_identifiers_from_neoantigens=set([n.patient_identifierforninself.neoantigens])patient_identifiers_from_patients=set([p.identifierforpinself.patients.values()])# check that there are no repeated neoantigensneoantigen_identifiers=[n.identifierforninself.neoantigens]iflen(neoantigen_identifiers)!=len(set(neoantigen_identifiers)):raiseNeofoxDataValidationException("There are repeated neoantigens!")# checks that no neoantigen is referring to an empty patientif""inpatient_identifiers_from_neoantigensorNoneinpatient_identifiers_from_neoantigens:raiseNeofoxDataValidationException("There are neoantigens missing a reference to a patient")# checks that there is no neoantigen referring to a non existing patientmissing_patient_identifiers=patient_identifiers_from_neoantigens.difference(patient_identifiers_from_patients)iflen(missing_patient_identifiers)>0:raiseNeofoxDataValidationException("There are neoantigens referring to missing patients: {}".format(missing_patient_identifiers))
[docs]defget_annotations(self)->List[NeoantigenAnnotations]:""" Loads epitope data (if file has been not imported to R; colnames need to be changed), adds data to class that are needed to calculate, calls epitope class --> determination of epitope properties, write to txt file """logger.info("Starting NeoFox annotations...")# feature calculation for each epitopefutures=[]start=time.time()forneoantigeninself.neoantigens:patient=self.patients.get(neoantigen.patient_identifier)logger.debug("Neoantigen: {}".format(neoantigen.to_json(indent=3)))logger.debug("Patient: {}".format(patient.to_json(indent=3)))futures.append(self.dask_client.submit(NeoFox.annotate_neoantigen,neoantigen,patient,self.reference_folder,self.configuration))annotations=self.dask_client.gather(futures)end=time.time()logger.info("Elapsed time for annotating {} neoantigens {} seconds".format(len(self.neoantigens),int(end-start)))returnannotations
[docs]@staticmethoddefannotate_neoantigen(neoantigen:Neoantigen,patient:Patient,reference_folder:ReferenceFolder,configuration:DependenciesConfiguration):logger.info("Starting neoantigen annotation: {}".format(neoantigen.identifier))start=time.time()annotation=NeoantigenAnnotator(reference_folder,configuration).get_annotation(neoantigen,patient)end=time.time()logger.info("Elapsed time for annotating neoantigen {}: {} seconds".format(neoantigen.identifier,int(end-start)))returnannotation
#!/usr/bin/env python## Copyright (c) 2020-2030 Translational Oncology at the Medical Center of the Johannes Gutenberg-University Mainz gGmbH.## This file is part of Neofox# (see https://github.com/tron-bioinformatics/neofox).## This program is free software: you can redistribute it and/or modify# it under the terms of the GNU General Public License as published by# the Free Software Foundation, either version 3 of the License, or# (at your option) any later version.## This program is distributed in the hope that it will be useful,# but WITHOUT ANY WARRANTY; without even the implied warranty of# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the# GNU General Public License for more details.## You should have received a copy of the GNU General Public License# along with this program. If not, see <http://www.gnu.org/licenses/>.#fromlogzeroimportloggerfromdatetimeimportdatetimeimportneofoxfromneofox.annotation_resources.uniprot.uniprotimportUniprotfromneofox.helpers.epitope_helperimportEpitopeHelperfromneofox.helpers.runnerimportRunnerfromneofox.MHC_predictors.MixMHCpred.mixmhc2predimportMixMhc2Predfromneofox.MHC_predictors.MixMHCpred.mixmhcpredimportMixMHCpredfromneofox.MHC_predictors.netmhcpan.combine_netmhcIIpan_pred_multiple_bindersimportBestAndMultipleBinderMhcIIfromneofox.MHC_predictors.netmhcpan.combine_netmhcpan_pred_multiple_bindersimportBestAndMultipleBinderfromneofox.published_features.differential_binding.amplitudeimportAmplitudefromneofox.published_features.differential_binding.differential_bindingimportDifferentialBindingfromneofox.published_features.Tcell_predictor.tcellpredictor_wrapperimportTcellPredictionfromneofox.published_features.dissimilarity_garnish.dissimilaritycalculatorimportDissimilarityCalculatorfromneofox.published_features.neoag.neoag_gbm_modelimportNeoagCalculatorfromneofox.published_features.neoantigen_fitness.neoantigen_fitnessimportNeoantigenFitnessCalculatorfromneofox.published_features.self_similarity.self_similarityimportSelfSimilarityCalculatorfromneofox.published_features.vaxrankimportvaxrankfromneofox.published_features.iedb_immunogenicity.iedbimportIEDBimmunogenicityfromneofox.published_features.expressionimportExpressionfromneofox.published_features.priority_scoreimportPriorityScorefromneofox.model.neoantigenimportPatient,Neoantigen,NeoantigenAnnotationsfromneofox.references.referencesimportReferenceFolder,DependenciesConfiguration
[docs]classNeoantigenAnnotator:def__init__(self,references:ReferenceFolder,configuration:DependenciesConfiguration):"""class to annotate neoantigens"""runner=Runner()self.dissimilarity_calculator=DissimilarityCalculator(runner=runner,configuration=configuration,proteome_db=references.proteome_db)self.neoantigen_fitness_calculator=NeoantigenFitnessCalculator(runner=runner,configuration=configuration,iedb=references.iedb)self.neoag_calculator=NeoagCalculator(runner=runner,configuration=configuration)self.netmhc2pan=BestAndMultipleBinderMhcII(runner=runner,configuration=configuration)self.mixmhc2=MixMhc2Pred(runner=runner,configuration=configuration)self.netmhcpan=BestAndMultipleBinder(runner=runner,configuration=configuration)self.mixmhc=MixMHCpred(runner=runner,configuration=configuration)self.available_alleles=references.get_available_alleles()self.uniprot=Uniprot(references.uniprot)self.tcell_predictor=TcellPrediction()self.self_similarity=SelfSimilarityCalculator()# NOTE: these resources do not read any file thus can be initialised fastself.differential_binding=DifferentialBinding()self.priority_score_calculator=PriorityScore()self.iedb_immunogenicity=IEDBimmunogenicity()self.differential_binding=DifferentialBinding()self.amplitude=Amplitude()
[docs]defget_annotation(self,neoantigen:Neoantigen,patient:Patient)->NeoantigenAnnotations:"""Calculate new epitope features and add to dictonary that stores all properties"""self._initialise_annotations(neoantigen)# decides which VAF to usevaf_rna=neoantigen.rna_variant_allele_frequencyifnotpatient.is_rna_available:logger.warning("Using the DNA VAF to estimate the RNA VAF as the patient does not have RNA available")# TODO: overwrite value in the neoantigen objectvaf_rna=neoantigen.dna_variant_allele_frequency# TODO: this is needed by the T cell predictor, move this construction inside by passing the neoantigensubstitution="{}{}{}".format(neoantigen.mutation.wild_type_aminoacid,neoantigen.mutation.position,neoantigen.mutation.mutated_aminoacid)# MHC binding independent featuresexpression_calculator=Expression(transcript_expression=neoantigen.rna_expression,vaf_rna=vaf_rna)self.annotations.annotations.extend(expression_calculator.get_annotations())sequence_not_in_uniprot=self.uniprot.is_sequence_not_in_uniprot(neoantigen.mutation.mutated_xmer)self.annotations.annotations.extend(self.uniprot.get_annotations(sequence_not_in_uniprot))# HLA I predictions: NetMHCpanself.netmhcpan.run(sequence_mut=neoantigen.mutation.mutated_xmer,sequence_wt=neoantigen.mutation.wild_type_xmer,mhc=patient.mhc1,available_mhc_alleles=self.available_alleles.get_available_mhc_i())self.annotations.annotations.extend(self.netmhcpan.get_annotations())# HLA II predictions: NetMHCIIpanself.netmhc2pan.run(sequence=neoantigen.mutation.mutated_xmer,sequence_reference=neoantigen.mutation.wild_type_xmer,mhc=patient.mhc2,available_mhc=self.available_alleles.get_available_mhc_ii())self.annotations.annotations.extend(self.netmhc2pan.get_annotations())# Amplitudeself.amplitude.run(netmhcpan=self.netmhcpan,netmhc2pan=self.netmhc2pan)self.annotations.annotations.extend(self.amplitude.get_annotations())# Neoantigen fitnessself.annotations.annotations.extend(self.neoantigen_fitness_calculator.get_annotations(self.netmhcpan,self.amplitude))# Differential Bindingself.annotations.annotations.extend(self.differential_binding.get_annotations_dai(self.netmhcpan))self.annotations.annotations.extend(self.differential_binding.get_annotations(self.netmhcpan,self.amplitude))self.annotations.annotations.extend(self.differential_binding.get_annotations_mhc2(self.netmhc2pan,self.amplitude))# T cell predictorself.annotations.annotations.extend(self.tcell_predictor.get_annotations(gene=neoantigen.transcript.gene,substitution=substitution,netmhcpan=self.netmhcpan))# self-similarityself.annotations.annotations.extend(self.self_similarity.get_annnotations(netmhcpan=self.netmhcpan))# number of mismatches and priority scoreself.annotations.annotations.extend(self.priority_score_calculator.get_annotations(netmhcpan=self.netmhcpan,vaf_transcr=vaf_rna,vaf_tum=neoantigen.dna_variant_allele_frequency,expr=neoantigen.rna_expression,mut_not_in_prot=sequence_not_in_uniprot))# neoag immunogenicity modelpeptide_variant_position=EpitopeHelper.position_of_mutation_epitope(wild_type=self.netmhcpan.best4_affinity_epitope_WT,mutation=self.netmhcpan.best4_affinity_epitope)self.annotations.annotations.append(self.neoag_calculator.get_annotation(sample_id=patient.identifier,netmhcpan=self.netmhcpan,peptide_variant_position=peptide_variant_position))# IEDB immunogenicityself.annotations.annotations.extend(self.iedb_immunogenicity.get_annotations(netmhcpan=self.netmhcpan,mhci_allele=self.netmhcpan.best4_affinity_allele))# MixMHCpredself.mixmhc.run(sequence_wt=neoantigen.mutation.wild_type_xmer,sequence_mut=neoantigen.mutation.mutated_xmer,mhc=patient.mhc1)self.annotations.annotations.extend(self.mixmhc.get_annotations())# MixMHC2predself.mixmhc2.run(mhc=patient.mhc2,sequence_wt=neoantigen.mutation.wild_type_xmer,sequence_mut=neoantigen.mutation.mutated_xmer)self.annotations.annotations.extend(self.mixmhc2.get_annotations())# dissimilarity to self-proteomeself.annotations.annotations.extend(self.dissimilarity_calculator.get_annotations(netmhcpan=self.netmhcpan))# vaxrankvaxrankscore=vaxrank.VaxRank()vaxrankscore.run(mutation_scores=self.netmhcpan.epitope_affinities,expression_score=expression_calculator.expression)self.annotations.annotations.extend(vaxrankscore.get_annotations())returnself.annotations
def_initialise_annotations(self,neoantigen):self.annotations=NeoantigenAnnotations()self.annotations.neoantigen_identifier=neoantigen.identifierself.annotations.annotator="Neofox"self.annotations.annotator_version=neofox.VERSIONself.annotations.timestamp="{:%Y%m%d%H%M%S%f}".format(datetime.now())# TODO: set the hash fro the resourcesself.annotations.annotations=[]
## Copyright (c) 2020-2030 Translational Oncology at the Medical Center of the Johannes Gutenberg-University Mainz gGmbH.## This file is part of Neofox# (see https://github.com/tron-bioinformatics/neofox).## This program is free software: you can redistribute it and/or modify# it under the terms of the GNU General Public License as published by# the Free Software Foundation, either version 3 of the License, or# (at your option) any later version.## This program is distributed in the hope that it will be useful,# but WITHOUT ANY WARRANTY; without even the implied warranty of# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the# GNU General Public License for more details.## You should have received a copy of the GNU General Public License# along with this program. If not, see <http://www.gnu.org/licenses/>.#
## Copyright (c) 2020-2030 Translational Oncology at the Medical Center of the Johannes Gutenberg-University Mainz gGmbH.## This file is part of Neofox# (see https://github.com/tron-bioinformatics/neofox).## This program is free software: you can redistribute it and/or modify# it under the terms of the GNU General Public License as published by# the Free Software Foundation, either version 3 of the License, or# (at your option) any later version.## This program is distributed in the hope that it will be useful,# but WITHOUT ANY WARRANTY; without even the implied warranty of# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the# GNU General Public License for more details.## You should have received a copy of the GNU General Public License# along with this program. If not, see <http://www.gnu.org/licenses/>.#importdataclassesimportenumimportinspectimportjsonimportstructimportsysfromabcimportABCfrombase64importb64encode,b64decodefromdatetimeimportdatetime,timedelta,timezonefromtypingimport(Any,AsyncGenerator,Callable,Collection,Dict,Generator,Iterable,List,Mapping,Optional,SupportsBytes,Tuple,Type,TypeVar,Union,get_type_hints,TYPE_CHECKING,)importgrpclib.constimportstringcasefrom.casingimportsafe_snake_caseifTYPE_CHECKING:fromgrpclib._protocolsimportIProtoMessagefromgrpclib.clientimportChannelfromgrpclib.metadataimportDeadlineifnot(sys.version_info.major==3andsys.version_info.minor>=7):# Apply backport of datetime.fromisoformat from 3.7frombackports.datetime_fromisoformatimportMonkeyPatchMonkeyPatch.patch_fromisoformat()# Proto 3 data typesTYPE_ENUM="enum"TYPE_BOOL="bool"TYPE_INT32="int32"TYPE_INT64="int64"TYPE_UINT32="uint32"TYPE_UINT64="uint64"TYPE_SINT32="sint32"TYPE_SINT64="sint64"TYPE_FLOAT="float"TYPE_DOUBLE="double"TYPE_FIXED32="fixed32"TYPE_SFIXED32="sfixed32"TYPE_FIXED64="fixed64"TYPE_SFIXED64="sfixed64"TYPE_STRING="string"TYPE_BYTES="bytes"TYPE_MESSAGE="message"TYPE_MAP="map"# Fields that use a fixed amount of space (4 or 8 bytes)FIXED_TYPES=[TYPE_FLOAT,TYPE_DOUBLE,TYPE_FIXED32,TYPE_SFIXED32,TYPE_FIXED64,TYPE_SFIXED64,]# Fields that are numerical 64-bit typesINT_64_TYPES=[TYPE_INT64,TYPE_UINT64,TYPE_SINT64,TYPE_FIXED64,TYPE_SFIXED64]# Fields that are efficiently packed whenPACKED_TYPES=[TYPE_ENUM,TYPE_BOOL,TYPE_INT32,TYPE_INT64,TYPE_UINT32,TYPE_UINT64,TYPE_SINT32,TYPE_SINT64,TYPE_FLOAT,TYPE_DOUBLE,TYPE_FIXED32,TYPE_SFIXED32,TYPE_FIXED64,TYPE_SFIXED64,]# Wire types# https://developers.google.com/protocol-buffers/docs/encoding#structureWIRE_VARINT=0WIRE_FIXED_64=1WIRE_LEN_DELIM=2WIRE_FIXED_32=5# Mappings of which Proto 3 types correspond to which wire types.WIRE_VARINT_TYPES=[TYPE_ENUM,TYPE_BOOL,TYPE_INT32,TYPE_INT64,TYPE_UINT32,TYPE_UINT64,TYPE_SINT32,TYPE_SINT64,]WIRE_FIXED_32_TYPES=[TYPE_FLOAT,TYPE_FIXED32,TYPE_SFIXED32]WIRE_FIXED_64_TYPES=[TYPE_DOUBLE,TYPE_FIXED64,TYPE_SFIXED64]WIRE_LEN_DELIM_TYPES=[TYPE_STRING,TYPE_BYTES,TYPE_MESSAGE,TYPE_MAP]# Protobuf datetimes start at the Unix Epoch in 1970 in UTC.defdatetime_default_gen():returndatetime(1970,1,1,tzinfo=timezone.utc)DATETIME_ZERO=datetime_default_gen()classCasing(enum.Enum):"""Casing constants for serialization."""CAMEL=stringcase.camelcaseSNAKE=stringcase.snakecaseclass_PLACEHOLDER:passPLACEHOLDER:Any=_PLACEHOLDER()@dataclasses.dataclass(frozen=True)classFieldMetadata:"""Stores internal metadata used for parsing & serialization."""# Protobuf field numbernumber:int# Protobuf type nameproto_type:str# Map information if the proto_type is a mapmap_types:Optional[Tuple[str,str]]=None# Groups several "one-of" fields togethergroup:Optional[str]=None# Describes the wrapped type (e.g. when using google.protobuf.BoolValue)wraps:Optional[str]=None@staticmethoddefget(field:dataclasses.Field)->"FieldMetadata":"""Returns the field metadata for a dataclass field."""returnfield.metadata["betterproto"]defdataclass_field(number:int,proto_type:str,*,map_types:Optional[Tuple[str,str]]=None,group:Optional[str]=None,wraps:Optional[str]=None,)->dataclasses.Field:"""Creates a dataclass field with attached protobuf metadata."""returndataclasses.field(default=PLACEHOLDER,metadata={"betterproto":FieldMetadata(number,proto_type,map_types,group,wraps)},)# Note: the fields below return `Any` to prevent type errors in the generated# data classes since the types won't match with `Field` and they get swapped# out at runtime. The generated dataclass variables are still typed correctly.defenum_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_ENUM,group=group)defbool_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_BOOL,group=group)defint32_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_INT32,group=group)defint64_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_INT64,group=group)defuint32_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_UINT32,group=group)defuint64_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_UINT64,group=group)defsint32_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_SINT32,group=group)defsint64_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_SINT64,group=group)deffloat_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_FLOAT,group=group)defdouble_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_DOUBLE,group=group)deffixed32_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_FIXED32,group=group)deffixed64_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_FIXED64,group=group)defsfixed32_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_SFIXED32,group=group)defsfixed64_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_SFIXED64,group=group)defstring_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_STRING,group=group)defbytes_field(number:int,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_BYTES,group=group)defmessage_field(number:int,group:Optional[str]=None,wraps:Optional[str]=None)->Any:returndataclass_field(number,TYPE_MESSAGE,group=group,wraps=wraps)defmap_field(number:int,key_type:str,value_type:str,group:Optional[str]=None)->Any:returndataclass_field(number,TYPE_MAP,map_types=(key_type,value_type),group=group)classEnum(int,enum.Enum):"""Protocol buffers enumeration base class. Acts like `enum.IntEnum`."""@classmethoddeffrom_string(cls,name:str)->int:"""Return the value which corresponds to the string name."""try:returncls.__members__[name]exceptKeyErrorase:raiseValueError(f"Unknown value {name} for enum {cls.__name__}")fromedef_pack_fmt(proto_type:str)->str:"""Returns a little-endian format string for reading/writing binary."""return{TYPE_DOUBLE:"<d",TYPE_FLOAT:"<f",TYPE_FIXED32:"<I",TYPE_FIXED64:"<Q",TYPE_SFIXED32:"<i",TYPE_SFIXED64:"<q",}[proto_type]defencode_varint(value:int)->bytes:"""Encodes a single varint value for serialization."""b:List[int]=[]ifvalue<0:value+=1<<64bits=value&0x7Fvalue>>=7whilevalue:b.append(0x80|bits)bits=value&0x7Fvalue>>=7returnbytes(b+[bits])def_preprocess_single(proto_type:str,wraps:str,value:Any)->bytes:"""Adjusts values before serialization."""ifproto_typein[TYPE_ENUM,TYPE_BOOL,TYPE_INT32,TYPE_INT64,TYPE_UINT32,TYPE_UINT64,]:returnencode_varint(value)elifproto_typein[TYPE_SINT32,TYPE_SINT64]:# Handle zig-zag encoding.ifvalue>=0:value=value<<1else:value=(value<<1)^(~0)returnencode_varint(value)elifproto_typeinFIXED_TYPES:returnstruct.pack(_pack_fmt(proto_type),value)elifproto_type==TYPE_STRING:returnvalue.encode("utf-8")elifproto_type==TYPE_MESSAGE:ifisinstance(value,datetime):# Convert the `datetime` to a timestamp message.seconds=int(value.timestamp())nanos=int(value.microsecond*1e3)value=_Timestamp(seconds=seconds,nanos=nanos)elifisinstance(value,timedelta):# Convert the `timedelta` to a duration message.total_ms=value//timedelta(microseconds=1)seconds=int(total_ms/1e6)nanos=int((total_ms%1e6)*1e3)value=_Duration(seconds=seconds,nanos=nanos)elifwraps:ifvalueisNone:returnb""value=_get_wrapper(wraps)(value=value)returnbytes(value)returnvaluedef_serialize_single(field_number:int,proto_type:str,value:Any,*,serialize_empty:bool=False,wraps:str="",)->bytes:"""Serializes a single field and value."""value=_preprocess_single(proto_type,wraps,value)output=b""ifproto_typeinWIRE_VARINT_TYPES:key=encode_varint(field_number<<3)output+=key+valueelifproto_typeinWIRE_FIXED_32_TYPES:key=encode_varint((field_number<<3)|5)output+=key+valueelifproto_typeinWIRE_FIXED_64_TYPES:key=encode_varint((field_number<<3)|1)output+=key+valueelifproto_typeinWIRE_LEN_DELIM_TYPES:iflen(value)orserialize_emptyorwraps:key=encode_varint((field_number<<3)|2)output+=key+encode_varint(len(value))+valueelse:raiseNotImplementedError(proto_type)returnoutputdefdecode_varint(buffer:bytes,pos:int,signed:bool=False)->Tuple[int,int]:""" Decode a single varint value from a byte buffer. Returns the value and the new position in the buffer. """result=0shift=0while1:b=buffer[pos]result|=(b&0x7F)<<shiftpos+=1ifnot(b&0x80):return(result,pos)shift+=7ifshift>=64:raiseValueError("Too many bytes when decoding varint.")@dataclasses.dataclass(frozen=True)classParsedField:number:intwire_type:intvalue:Anyraw:bytesdefparse_fields(value:bytes)->Generator[ParsedField,None,None]:i=0whilei<len(value):start=inum_wire,i=decode_varint(value,i)number=num_wire>>3wire_type=num_wire&0x7decoded:Any=Noneifwire_type==0:decoded,i=decode_varint(value,i)elifwire_type==1:decoded,i=value[i:i+8],i+8elifwire_type==2:length,i=decode_varint(value,i)decoded=value[i:i+length]i+=lengthelifwire_type==5:decoded,i=value[i:i+4],i+4yieldParsedField(number=number,wire_type=wire_type,value=decoded,raw=value[start:i])# Bound type variable to allow methods to return `self` of subclassesT=TypeVar("T",bound="Message")classProtoClassMetadata:cls:Type["Message"]def__init__(self,cls:Type["Message"]):self.cls=clsby_field={}by_group={}forfieldindataclasses.fields(cls):meta=FieldMetadata.get(field)ifmeta.group:# This is part of a one-of group.by_field[field.name]=meta.groupby_group.setdefault(meta.group,set()).add(field)self.oneof_group_by_field=by_fieldself.oneof_field_by_group=by_groupself.init_default_gen()self.init_cls_by_field()definit_default_gen(self):default_gen={}forfieldindataclasses.fields(self.cls):meta=FieldMetadata.get(field)default_gen[field.name]=self.cls._get_field_default_gen(field,meta)self.default_gen=default_gendefinit_cls_by_field(self):field_cls={}forfieldindataclasses.fields(self.cls):meta=FieldMetadata.get(field)ifmeta.proto_type==TYPE_MAP:assertmeta.map_typeskt=self.cls._cls_for(field,index=0)vt=self.cls._cls_for(field,index=1)Entry=dataclasses.make_dataclass("Entry",[("key",kt,dataclass_field(1,meta.map_types[0])),("value",vt,dataclass_field(2,meta.map_types[1])),],bases=(Message,),)field_cls[field.name]=Entryfield_cls[field.name+".value"]=vtelse:field_cls[field.name]=self.cls._cls_for(field)self.cls_by_field=field_clsclassMessage(ABC):""" A protobuf message base class. Generated code will inherit from this and register the message fields which get used by the serializers and parsers to go between Python, binary and JSON protobuf message representations. """_serialized_on_wire:bool_unknown_fields:bytes_group_map:Dict[str,dict]def__post_init__(self)->None:# Keep track of whether every field was defaultall_sentinel=True# Set a default value for each field in the class after `__init__` has# already been run.group_map:Dict[str,dataclasses.Field]={}forfieldindataclasses.fields(self):meta=FieldMetadata.get(field)ifmeta.group:group_map.setdefault(meta.group)ifgetattr(self,field.name)!=PLACEHOLDER:# Skip anything not set to the sentinel valueall_sentinel=Falseifmeta.group:# This was set, so make it the selected value of the one-of.group_map[meta.group]=fieldcontinuesetattr(self,field.name,self._get_field_default(field,meta))# Now that all the defaults are set, reset it!self.__dict__["_serialized_on_wire"]=notall_sentinelself.__dict__["_unknown_fields"]=b""self.__dict__["_group_map"]=group_mapdef__setattr__(self,attr:str,value:Any)->None:ifattr!="_serialized_on_wire":# Track when a field has been set.self.__dict__["_serialized_on_wire"]=Trueifhasattr(self,"_group_map"):# __post_init__ had already runifattrinself._betterproto.oneof_group_by_field:group=self._betterproto.oneof_group_by_field[attr]forfieldinself._betterproto.oneof_field_by_group[group]:iffield.name==attr:self._group_map[group]=fieldelse:super().__setattr__(field.name,self._get_field_default(field,FieldMetadata.get(field)),)super().__setattr__(attr,value)@propertydef_betterproto(self):""" Lazy initialize metadata for each protobuf class. It may be initialized multiple times in a multi-threaded environment, but that won't affect the correctness. """meta=getattr(self.__class__,"_betterproto_meta",None)ifnotmeta:meta=ProtoClassMetadata(self.__class__)self.__class__._betterproto_meta=metareturnmetadef__bytes__(self)->bytes:""" Get the binary encoded Protobuf representation of this instance. """output=b""forfieldindataclasses.fields(self):meta=FieldMetadata.get(field)value=getattr(self,field.name)ifvalueisNone:# Optional items should be skipped. This is used for the Google# wrapper types.continue# Being selected in a a group means this field is the one that is# currently set in a `oneof` group, so it must be serialized even# if the value is the default zero value.selected_in_group=Falseifmeta.groupandself._group_map[meta.group]==field:selected_in_group=Trueserialize_empty=Falseifisinstance(value,Message)andvalue._serialized_on_wire:# Empty messages can still be sent on the wire if they were# set (or received empty).serialize_empty=Trueifvalue==self._get_field_default(field,meta)andnot(selected_in_grouporserialize_empty):# Default (zero) values are not serialized. Two exceptions are# if this is the selected oneof item or if we know we have to# serialize an empty message (i.e. zero value was explicitly# set by the user).continueifisinstance(value,list):ifmeta.proto_typeinPACKED_TYPES:# Packed lists look like a length-delimited field. First,# preprocess/encode each value into a buffer and then# treat it like a field of raw bytes.buf=b""foriteminvalue:buf+=_preprocess_single(meta.proto_type,"",item)output+=_serialize_single(meta.number,TYPE_BYTES,buf)else:foriteminvalue:output+=_serialize_single(meta.number,meta.proto_type,item,wraps=meta.wrapsor"")elifisinstance(value,dict):fork,vinvalue.items():assertmeta.map_typessk=_serialize_single(1,meta.map_types[0],k)sv=_serialize_single(2,meta.map_types[1],v)output+=_serialize_single(meta.number,meta.proto_type,sk+sv)else:output+=_serialize_single(meta.number,meta.proto_type,value,serialize_empty=serialize_empty,wraps=meta.wrapsor"",)returnoutput+self._unknown_fields# For compatibility with other librariesSerializeToString=__bytes__@classmethoddef_type_hint(cls,field_name:str)->Type:module=inspect.getmodule(cls)type_hints=get_type_hints(cls,vars(module))returntype_hints[field_name]@classmethoddef_cls_for(cls,field:dataclasses.Field,index:int=0)->Type:"""Get the message class for a field from the type hints."""field_cls=cls._type_hint(field.name)ifhasattr(field_cls,"__args__")andindex>=0:field_cls=field_cls.__args__[index]returnfield_clsdef_get_field_default(self,field:dataclasses.Field,meta:FieldMetadata)->Any:returnself._betterproto.default_gen[field.name]()@classmethoddef_get_field_default_gen(cls,field:dataclasses.Field,meta:FieldMetadata)->Any:t=cls._type_hint(field.name)ifhasattr(t,"__origin__"):ift.__origin__in(dict,Dict):# This is some kind of map (dict in Python).returndictelift.__origin__in(list,List):# This is some kind of list (repeated) field.returnlistelift.__origin__==Unionandt.__args__[1]==type(None):# This is an optional (wrapped) field. For setting the default we# really don't care what kind of field it is.returntype(None)else:returntelifissubclass(t,Enum):# Enums always default to zero.returnintelift==datetime:# Offsets are relative to 1970-01-01T00:00:00Zreturndatetime_default_genelse:# This is either a primitive scalar or another message type. Calling# it should result in its zero value.returntdef_postprocess_single(self,wire_type:int,meta:FieldMetadata,field:dataclasses.Field,value:Any)->Any:"""Adjusts values after parsing."""ifwire_type==WIRE_VARINT:ifmeta.proto_typein[TYPE_INT32,TYPE_INT64]:bits=int(meta.proto_type[3:])value=value&((1<<bits)-1)signbit=1<<(bits-1)value=int((value^signbit)-signbit)elifmeta.proto_typein[TYPE_SINT32,TYPE_SINT64]:# Undo zig-zag encodingvalue=(value>>1)^(-(value&1))elifmeta.proto_type==TYPE_BOOL:# Booleans use a varint encoding, so convert it to true/false.value=value>0elifwire_typein[WIRE_FIXED_32,WIRE_FIXED_64]:fmt=_pack_fmt(meta.proto_type)value=struct.unpack(fmt,value)[0]elifwire_type==WIRE_LEN_DELIM:ifmeta.proto_type==TYPE_STRING:value=value.decode("utf-8")elifmeta.proto_type==TYPE_MESSAGE:cls=self._betterproto.cls_by_field[field.name]ifcls==datetime:value=_Timestamp().parse(value).to_datetime()elifcls==timedelta:value=_Duration().parse(value).to_timedelta()elifmeta.wraps:# This is a Google wrapper value message around a single# scalar type.value=_get_wrapper(meta.wraps)().parse(value).valueelse:value=cls().parse(value)value._serialized_on_wire=Trueelifmeta.proto_type==TYPE_MAP:value=self._betterproto.cls_by_field[field.name]().parse(value)returnvaluedefparse(self:T,data:bytes)->T:""" Parse the binary encoded Protobuf into this message instance. This returns the instance itself and is therefore assignable and chainable. """fields={f.metadata["betterproto"].number:fforfindataclasses.fields(self)}forparsedinparse_fields(data):ifparsed.numberinfields:field=fields[parsed.number]meta=FieldMetadata.get(field)value:Anyif(parsed.wire_type==WIRE_LEN_DELIMandmeta.proto_typeinPACKED_TYPES):# This is a packed repeated field.pos=0value=[]whilepos<len(parsed.value):ifmeta.proto_typein["float","fixed32","sfixed32"]:decoded,pos=parsed.value[pos:pos+4],pos+4wire_type=WIRE_FIXED_32elifmeta.proto_typein["double","fixed64","sfixed64"]:decoded,pos=parsed.value[pos:pos+8],pos+8wire_type=WIRE_FIXED_64else:decoded,pos=decode_varint(parsed.value,pos)wire_type=WIRE_VARINTdecoded=self._postprocess_single(wire_type,meta,field,decoded)value.append(decoded)else:value=self._postprocess_single(parsed.wire_type,meta,field,parsed.value)current=getattr(self,field.name)ifmeta.proto_type==TYPE_MAP:# Value represents a single key/value pair entry in the map.current[value.key]=value.valueelifisinstance(current,list)andnotisinstance(value,list):current.append(value)else:setattr(self,field.name,value)else:self._unknown_fields+=parsed.rawreturnself# For compatibility with other libraries.@classmethoddefFromString(cls:Type[T],data:bytes)->T:returncls().parse(data)defto_dict(self,casing:Casing=Casing.CAMEL,include_default_values:bool=False)->dict:""" Returns a dict representation of this message instance which can be used to serialize to e.g. JSON. Defaults to camel casing for compatibility but can be set to other modes. `include_default_values` can be set to `True` to include default values of fields. E.g. an `int32` type field with `0` value will not be in returned dict if `include_default_values` is set to `False`. """output:Dict[str,Any]={}forfieldindataclasses.fields(self):meta=FieldMetadata.get(field)v=getattr(self,field.name)cased_name=casing(field.name).rstrip("_")# type: ignoreifmeta.proto_type=="message":ifisinstance(v,datetime):ifv!=DATETIME_ZEROorinclude_default_values:output[cased_name]=_Timestamp.timestamp_to_json(v)elifisinstance(v,timedelta):ifv!=timedelta(0)orinclude_default_values:output[cased_name]=_Duration.delta_to_json(v)elifmeta.wraps:ifvisnotNoneorinclude_default_values:output[cased_name]=velifisinstance(v,list):# Convert each item.v=[i.to_dict(casing,include_default_values)foriinv]ifvorinclude_default_values:output[cased_name]=velse:ifv._serialized_on_wireorinclude_default_values:output[cased_name]=v.to_dict(casing,include_default_values)elifmeta.proto_type=="map":forkinv:ifhasattr(v[k],"to_dict"):v[k]=v[k].to_dict(casing,include_default_values)ifvorinclude_default_values:output[cased_name]=velifv!=self._get_field_default(field,meta)orinclude_default_values:ifmeta.proto_typeinINT_64_TYPES:ifisinstance(v,list):output[cased_name]=[str(n)forninv]else:output[cased_name]=str(v)elifmeta.proto_type==TYPE_BYTES:ifisinstance(v,list):output[cased_name]=[b64encode(b).decode("utf8")forbinv]else:output[cased_name]=b64encode(v).decode("utf8")elifmeta.proto_type==TYPE_ENUM:enum_values=list(self._betterproto.cls_by_field[field.name])# type: ignoreifisinstance(v,list):output[cased_name]=[enum_values[e].nameforeinv]else:output[cased_name]=enum_values[v].nameelse:output[cased_name]=vreturnoutputdeffrom_dict(self:T,value:dict)->T:""" Parse the key/value pairs in `value` into this message instance. This returns the instance itself and is therefore assignable and chainable. """self._serialized_on_wire=Truefields_by_name={f.name:fforfindataclasses.fields(self)}forkeyinvalue:snake_cased=safe_snake_case(key)ifsnake_casedinfields_by_name:field=fields_by_name[snake_cased]meta=FieldMetadata.get(field)ifvalue[key]isnotNone:ifmeta.proto_type=="message":v=getattr(self,field.name)ifisinstance(v,list):cls=self._betterproto.cls_by_field[field.name]foriinrange(len(value[key])):v.append(cls().from_dict(value[key][i]))elifisinstance(v,datetime):v=datetime.fromisoformat(value[key].replace("Z","+00:00"))setattr(self,field.name,v)elifisinstance(v,timedelta):v=timedelta(seconds=float(value[key][:-1]))setattr(self,field.name,v)elifmeta.wraps:setattr(self,field.name,value[key])else:v.from_dict(value[key])elifmeta.map_typesandmeta.map_types[1]==TYPE_MESSAGE:v=getattr(self,field.name)cls=self._betterproto.cls_by_field[field.name+".value"]forkinvalue[key]:v[k]=cls().from_dict(value[key][k])else:v=value[key]ifmeta.proto_typeinINT_64_TYPES:ifisinstance(value[key],list):v=[int(n)forninvalue[key]]else:v=int(value[key])elifmeta.proto_type==TYPE_BYTES:ifisinstance(value[key],list):v=[b64decode(n)forninvalue[key]]else:v=b64decode(value[key])elifmeta.proto_type==TYPE_ENUM:enum_cls=self._betterproto.cls_by_field[field.name]ifisinstance(v,list):v=[enum_cls.from_string(e)foreinv]elifisinstance(v,str):v=enum_cls.from_string(v)ifvisnotNone:setattr(self,field.name,v)returnselfdefto_json(self,indent:Union[None,int,str]=None)->str:"""Returns the encoded JSON representation of this message instance."""returnjson.dumps(self.to_dict(),indent=indent)deffrom_json(self:T,value:Union[str,bytes])->T:""" Parse the key/value pairs in `value` into this message instance. This returns the instance itself and is therefore assignable and chainable. """returnself.from_dict(json.loads(value))defserialized_on_wire(message:Message)->bool:""" True if this message was or should be serialized on the wire. This can be used to detect presence (e.g. optional wrapper message) and is used internally during parsing/serialization. """returnmessage._serialized_on_wiredefwhich_one_of(message:Message,group_name:str)->Tuple[str,Any]:"""Return the name and value of a message's one-of field group."""field=message._group_map.get(group_name)ifnotfield:return("",None)return(field.name,getattr(message,field.name))@dataclasses.dataclassclass_Duration(Message):# Signed seconds of the span of time. Must be from -315,576,000,000 to# +315,576,000,000 inclusive. Note: these bounds are computed from: 60# sec/min * 60 min/hr * 24 hr/day * 365.25 days/year * 10000 yearsseconds:int=int64_field(1)# Signed fractions of a second at nanosecond resolution of the span of time.# Durations less than one second are represented with a 0 `seconds` field and# a positive or negative `nanos` field. For durations of one second or more,# a non-zero value for the `nanos` field must be of the same sign as the# `seconds` field. Must be from -999,999,999 to +999,999,999 inclusive.nanos:int=int32_field(2)defto_timedelta(self)->timedelta:returntimedelta(seconds=self.seconds,microseconds=self.nanos/1e3)@staticmethoddefdelta_to_json(delta:timedelta)->str:parts=str(delta.total_seconds()).split(".")iflen(parts)>1:whilelen(parts[1])notin[3,6,9]:parts[1]=parts[1]+"0"return".".join(parts)+"s"@dataclasses.dataclassclass_Timestamp(Message):# Represents seconds of UTC time since Unix epoch 1970-01-01T00:00:00Z. Must# be from 0001-01-01T00:00:00Z to 9999-12-31T23:59:59Z inclusive.seconds:int=int64_field(1)# Non-negative fractions of a second at nanosecond resolution. Negative# second values with fractions must still have non-negative nanos values that# count forward in time. Must be from 0 to 999,999,999 inclusive.nanos:int=int32_field(2)defto_datetime(self)->datetime:ts=self.seconds+(self.nanos/1e9)returndatetime.fromtimestamp(ts,tz=timezone.utc)@staticmethoddeftimestamp_to_json(dt:datetime)->str:nanos=dt.microsecond*1e3copy=dt.replace(microsecond=0,tzinfo=None)result=copy.isoformat()if(nanos%1e9)==0:# If there are 0 fractional digits, the fractional# point '.' should be omitted when serializing.returnresult+"Z"if(nanos%1e6)==0:# Serialize 3 fractional digits.returnresult+".%03dZ"%(nanos/1e6)if(nanos%1e3)==0:# Serialize 6 fractional digits.returnresult+".%06dZ"%(nanos/1e3)# Serialize 9 fractional digits.returnresult+".%09dZ"%nanosclass_WrappedMessage(Message):""" Google protobuf wrapper types base class. JSON representation is just the value itself. """value:Anydefto_dict(self,casing:Casing=Casing.CAMEL)->Any:returnself.valuedeffrom_dict(self:T,value:Any)->T:ifvalueisnotNone:self.value=valuereturnself@dataclasses.dataclassclass_BoolValue(_WrappedMessage):value:bool=bool_field(1)@dataclasses.dataclassclass_Int32Value(_WrappedMessage):value:int=int32_field(1)@dataclasses.dataclassclass_UInt32Value(_WrappedMessage):value:int=uint32_field(1)@dataclasses.dataclassclass_Int64Value(_WrappedMessage):value:int=int64_field(1)@dataclasses.dataclassclass_UInt64Value(_WrappedMessage):value:int=uint64_field(1)@dataclasses.dataclassclass_FloatValue(_WrappedMessage):value:float=float_field(1)@dataclasses.dataclassclass_DoubleValue(_WrappedMessage):value:float=double_field(1)@dataclasses.dataclassclass_StringValue(_WrappedMessage):value:str=string_field(1)@dataclasses.dataclassclass_BytesValue(_WrappedMessage):value:bytes=bytes_field(1)def_get_wrapper(proto_type:str)->Type:"""Get the wrapper message class for a wrapped type."""return{TYPE_BOOL:_BoolValue,TYPE_INT32:_Int32Value,TYPE_UINT32:_UInt32Value,TYPE_INT64:_Int64Value,TYPE_UINT64:_UInt64Value,TYPE_FLOAT:_FloatValue,TYPE_DOUBLE:_DoubleValue,TYPE_STRING:_StringValue,TYPE_BYTES:_BytesValue,}[proto_type]_Value=Union[str,bytes]_MetadataLike=Union[Mapping[str,_Value],Collection[Tuple[str,_Value]]]classServiceStub(ABC):""" Base class for async gRPC service stubs. """def__init__(self,channel:"Channel",*,timeout:Optional[float]=None,deadline:Optional["Deadline"]=None,metadata:Optional[_MetadataLike]=None,)->None:self.channel=channelself.timeout=timeoutself.deadline=deadlineself.metadata=metadatadef__resolve_request_kwargs(self,timeout:Optional[float],deadline:Optional["Deadline"],metadata:Optional[_MetadataLike],):return{"timeout":self.timeoutiftimeoutisNoneelsetimeout,"deadline":self.deadlineifdeadlineisNoneelsedeadline,"metadata":self.metadataifmetadataisNoneelsemetadata,}asyncdef_unary_unary(self,route:str,request:"IProtoMessage",response_type:Type[T],*,timeout:Optional[float]=None,deadline:Optional["Deadline"]=None,metadata:Optional[_MetadataLike]=None,)->T:"""Make a unary request and return the response."""asyncwithself.channel.request(route,grpclib.const.Cardinality.UNARY_UNARY,type(request),response_type,**self.__resolve_request_kwargs(timeout,deadline,metadata),)asstream:awaitstream.send_message(request,end=True)response=awaitstream.recv_message()assertresponseisnotNonereturnresponseasyncdef_unary_stream(self,route:str,request:"IProtoMessage",response_type:Type[T],*,timeout:Optional[float]=None,deadline:Optional["Deadline"]=None,metadata:Optional[_MetadataLike]=None,)->AsyncGenerator[T,None]:"""Make a unary request and return the stream response iterator."""asyncwithself.channel.request(route,grpclib.const.Cardinality.UNARY_STREAM,type(request),response_type,**self.__resolve_request_kwargs(timeout,deadline,metadata),)asstream:awaitstream.send_message(request,end=True)asyncformessageinstream:yieldmessage
class neofox.annotator.NeoantigenAnnotator(references:neofox.references.references.ReferenceFolder, configuration:neofox.references.references.DependenciesConfiguration)[source]¶
Loads epitope data (if file has been not imported to R; colnames need to be changed), adds data to class that are needed to calculate,
calls epitope class –> determination of epitope properties,
write to txt file
class neofox.model.neoantigen.Mhc2GeneName(value)[source]¶
Bases: betterproto.Enum
*Valid names for MHC II classic genes.DRA is not included in this list as
it does not have much variability in the population and for our purpose
isconsidered constant.
*MHC allele representation. It does not include non synonymous changes to
the sequence, changes in the non coding regionor changes in expression. See
http://hla.alleles.org/nomenclature/naming.html for details
class neofox.model.neoantigen.Mutation(position:int=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>, wild_type_xmer:str=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>, wild_type_aminoacid:str=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>, mutated_xmer:str=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>, mutated_aminoacid:str=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>, left_flanking_region:str=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>, size_left_flanking_region:int=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>, right_flanking_region:str=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>, size_right_flanking_region:int=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>)[source]¶
class neofox.model.neoantigen.Transcript(identifier:str=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>, assembly:str=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>, gene:str=<betterproto._PLACEHOLDER object at 0x7f55e71c67b8>)[source]¶
transcript.assembly - the assembly of the reference genome (only hg19 is supported)
transcript.gene - the HGMC gene symbol
transcript.identifier - the UCSC transcript identifier including the version number
mutation.leftFlankingRegion - the aminoacids flanking the mutation on the left (in IUPAC one letter symbols)
mutation.mutatedAminoacid - the mutated aminoacid (IUPAC 1 or 3 letters respecting casing, eg: A and Ala)
mutation.position - the 1 based position of the mutation in the protein
mutation.rightFlankingRegion - the aminoacids flanking the mutation on the left (in IUPAC one letter symbols)
mutation.wildTypeAminoacid - the wild type aminoacid (IUPAC 1 or 3 letters respecting casing, eg: A and Ala)
patientIdentifier - the patient identifier
rnaExpression - the transcript RNA expression (optional)
rnaVariantAlleleFrequency - the variant allele frequency calculated from the RNA (optional, this will be estimated using the dnaVariantAlleleFrequency if not available)
dnaVariantAlleleFrequency - the variant allele frequency calculated from the DNA (optional)
The output data is returned in a short wide tab separated values file (--with-short-wide-table). Optionally, it can be provided in a tall skinny tab separated values file (--with-tall-skinny-table) or in JSON (--with-json).
